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Abstract. - We study step bunching under conditions of attachment/detachment limited 
kinetics in the presence of a deposition or sublimation flux, which leads to bunch motion. 
Analysis of the discrete step dynamics reveals that the bunch velocity is inversely proportional 
to the bunch size for general step-step interactions. The shape of steadily moving bunches is 
studied within a continuum theory, and analytic expressions for the bunch profile are derived. 
Scaling laws obtained previously for non-moving bunches are recovered asymptotically, but 
singularities of the static theory are removed and strong corrections to scaling are found. The 
size of the largest terrace between two bunches is identified as a central scaling parameter. Our 
theory applies to a large class of bunching instabilities, including sublimation with attachment 
asymmetry and surface electromigration in the presence of sublimation or growth. 



Introduction. - There is much current interest in exploiting morphological instabilities 
to form periodic nanoscale patterns at crystal surfaces [1,2]. Because of their natural in-plane 
anisotropy, vicinal surfaces [3] prepared at a miscut relative to a high symmetry orientation 
provide ideal substrates for the formation of ripple patterns parallel or perpendicular to the 
mean step orientation [4, 5] . Here we specifically consider patterns formed by step bunching, 
the process in which a train of initially equally spaced (straight) steps splits into regions of 
densely packed steps (step bunches), and almost flat regions [6,7]. Bunched semiconductor 
surfaces are promising templates for the growth of metallic nanowires [8] . 

Step bunching can be induced by growth [9,10], sublimation [11], or surface migration of 
adatoms driven by an electric current [4, 12-18]. The common feature of the different insta- 
bility mechanisms [7] is that they break the symmetry between the ascending and descending 
steps bordering a terrace. The appearance of step bunches thus provides information about 
the asymmetry of the attachment/detachment processes at the steps, as well as about the 
direction of current-induced adatom migration. Once formed, the shape of a bunch is deter- 
mined by the balance between the destabilizing forces and the repulsive step-step interactions 
that act to regularize the step train. As a result, the bunch shape displays characteristic 
scaling laws relating e.g. its slope and width to the number of steps in the bunch [13,14]. 
These scaling laws are used in the interpretation of experiments to extract the functional form 
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of the step interactions as well as material parameters such as the step interaction strength 
and the electromigration force [17, 18]. 

The large scale properties of step bunches are captured by continuum evolution equations 
for the surface profile [6], which can be derived from the underlying discrete step dynamics 
in a systematic manner [11]. The analysis of static (time-independent) solutions of these 
equations leads to scaling laws which are in reasonable agreement with numerical simulation 
of the discrete step dynamics [11]. However, in the presence of a non- vanishing sublimation 
or growth flux, step bunches are moving objects. Because of the high temperatures involved, 
sublimation - and hence, bunch motion - is significant also in electromigration experiments, 
where it is not the primary cause of bunching [12, 15]. 

In this Letter we show that bunch motion alters the shape and scaling properties of bunches 
in a fundamental way. It removes the artificial symmetry between the in-flow and out-flow 
regions (in which steps move into and out of the bunch, respectively) and the concomitant 
singularities of the static solutions at the bunch edges [11]. We show that the lateral speed 
of a bunch is inversely proportional to its height for a large class of models, and we identify 
the size of the largest terrace l mayi as a natural scaling parameter, in terms of which other 
important bunch characteristics are expressed in a simple way. The maximal terrace size 
^max is uniquely defined, in contrast to the number of steps in the bunch, which requires a 
convention to decide which steps belong to it, and it is directly accessible experimentally by 
means of reflection electron microscopy (REM) [16]. 

Discrete model. We consider a system of non-transparent steps [7] described on the 
discrete level by the equations of motion 

= — ( x i+l - Xi) + ( x i ~ x i-l) + U (2/i - fi-i - f i+ i) (1) 

for the step positions Xi(t), where the time scale has been normalized to the growth or sub- 
limation flux. The parameter b governs the asymmetry between ascending and descending 
steps, relative to the mean step velocity. The linear form of the first two terms on the right 
hand side of JIJ is characteristic of slow attachment/detachment kinetics, and applies equally 
to step bunching induced by sublimation, growth or surface electromigration [11, 14]; here we 
will assume a sublimating step train going uphill in the +x direction. The last term on the 
right hand side of (IJ represents stabilizing step-step interactions of strength U. In the usual 
case of entropic or dipolar elastic interactions 

,,_k_L_r , -k_L^ 

\Xi Xi—i/ \Xi-\-i X' L 

where v — 2 and I is the average terrace length [3] . Explicit expressions for b and U in terms 
of physical parameters are given below in (|15(l . 

For b > 0, JU leads to an instability of the equally spaced step configuration x% — Xi-\ = I 
and its segregation into step bunches separated by flat regions. The bunches coarsen slowly in 
time by coalescence. We are interested in the final regime of coarsening with a few big bunches 
left in the system. In this regime, one can study a periodic array of identical bunches, each 
containing M steps, which satisfy (JTJ with i = 1,2, ...M and the helicoidal boundary conditions 
Xm+i = xi+Ml. It is convenient to consider the comoving step coordinates yi{i) = Xi(t)—lt, in 
which the center of mass of the step configuation does not move. In this frame, the stationary 
trajectory of a step is a periodic function with some (unknown) period r, yi{t) = yi(t + r). 
Stationarity implies that every step follows the same trajectory, up to a space and time shift, 
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according to yi+ s (t) = yi(t + ts/M) + Is, with s = 1, 2, ...M — 1. Inserting this into and 
setting A(i) = — Ui{t) we obtain an equation for the stationary step trajectory (in the 

following we omit subscripts) 



dt 



i 



-A(t) + ^A(t-^)-l + U 



2 /w -/<*-£) 



(3) 



This is a differential-difference equation for two periodic functions y(t) and /(t), which for 
the time being will be treated as independent. Expanding the functions in Fourier series with 
frequencies w n = 2-kti/t and coefficients Y n and F ni respectively, we obtain from |J3j 



sin ■ 



2im 
W 



27m \ 

■IV ( I -cos — ) Im[F„/y n 



(4) 



Since J3J is valid for any n, we can choose n/M <C 1 for large M and expand (J2J to obtain the 
expression r = M — 2ttU lm[nF n /Y n ] + 0(M ), which in fact determines the dependence of 
the bunch velocity on M: In the laboratory frame we have Xi(t + r/M) = Xi + i (t) + Z(r/M — 1), 
which implies that the whole step configuration shifts by l{r/M — I) to the right in time r/M . 
Hence the lateral bunch speed v is given by 



j/l=l- M/t = k/M + o(M" 1 ), 



(5) 



where k = — 2TrUlm[nF n /Y n ] has to be determined self-consistently for a given form of step- 
step interaction. When / has the usual form (J2J with v — 2, we find numerically that K is 
a constant proportional to the asymmetry, k ps 36, provided that the asymmetry is not too 
large, b < 0.5. 

In the following we will see that even though the bunch velocity decreases with increasing 
bunch size, it cannot be neglected. According to a scaling argument due to Chernov [9,13,19], 
the scaling v ~ M^ 1 implies that the average bunch size should increase with time as 
which is consistent with experiments [12,18] and numerous discrete simulation [14,15]. 



Continuum theory. 
dynamics reads [11] 



The continuum evolution equation corresponding to the discrete 



dh d 
dt dx 



M 

2m 6m 3 dx 2m dx 2 



"o 



dm 3Ul 3 d 2 (m 2 ) 



fro = 0, 



(6) 



where h(x,t) is the surface profile, m(x, t) — dh/dx > is the slope, and ho denotes the 
height of a single step. A periodic array of bunches moving at lateral speed v is obtained by 
the travelling wave ansatz [9] h(x,t) — h(£) + fit — hot, where £ = x — vt and the function 
h(£) satisfies the boundary condition + Ml) = h(£) + Mho- Inserting this into © we find 
that the vertical excess speed Q is related to v by Q — vmo, where mo = ho /I is the average 
slope. Integrating once the ordinary differential equation for then becomes 



n(t + z -h) + 



i 



l 

m 



ZU 



6m 3 2m 



2y/ 



(m 2 ) 



0. 



(7) 



We denote derivatives by primes, and measure lateral distances in units of I and heights 
in units of ho- The phase shift £o is a constant of integration satisfying the condition 
J Q M m(£) — £ — £o] d£ — 0, which follows by multiplying by m, integrating over the 
bunch period, and using the boundary conditions. For future reference we note that, for small 
b, tt sa 3b/ M because of ©. 
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Fig. 2 



Fig. 1 - Shape of a moving bunch h(£) computed from the discrete dynamics (J^J [symbols], and the 
continuous evolution equation Q [full line], for 128 steps with b = 3/17 ~ 0.17647, and U = 0.569. 
There are no fitting parameters. 

Fig. 2 - Slope profile derived from Fig0 for a moving bunch (full curve) and a static bunch (bold 
dashed curve). The slope of the moving bunch increases abruptly in the inflow region and decreases 
gradually in the outflow region. 



In Fig^ a numerical solution of Q obtained via a shooting method is compared to the 
discrete step dynamics, showing excellent agreement. The continuum description is generally 
found to work very well, provided the asymmetry 6 is sufficiently small. In Fig|21we compare 
the corresponding slope profile to the time-independent solution derived in [11] by setting 
the terms inside the square brackets in JSJ to zero and neglecting the symmetry-breaking 
term (h^/Qm^dm/ dx. As was anticipated in [11], the moving bunch is distinctly asym- 
metric. Moreover the moving solution extends smoothly over the whole x-axis, whereas the 
static solution has finite support due to singularities at the bunch edges. In the following 
we find analytically the asymptotics of the inflow and outflow regions of the moving bunch, 
corresponding to the extreme left and right parts of FigQ](see also Fig|2J). 

Outflow region: m(£) <C 1, m'(^) < 0. - In this region the steps are far apart and their 
interaction U is negligible. It is convenient to perform a Lagrange transform from the function 
h(£) to its inverse £ = £(h). After some algebra Q with U = then reduces to the linear 
equation 

f " - 36£' + 6Qf = 6Q(h - f„) - 36, (8) 
which has the general solution 

t = h-to+Ci exp (Ai/i) + C 2 cxp (X 2 h) (9) 

with Ai,2 = (36/2) [1 ± y/l — 8S7/36 2 ]. Fixing the boundary conditions so that the point of 
minimal slope min^ m(£) = e is located at h = 0, we have two boundary conditions £(0) = 
and £'(0) = l/e to determine C\,C 2 . Recalling that fl w 36/M, we see that Ai — > 36 
and A2 — * for large bunches, so that @ becomes a pure exponential, corresponding to a 
slope profile m(£) w l/(36£). Similar behavior was found in a model with short-range step 
interaction, however in that case SI — * 36 2 /8 and Ai — > A2 for large bunches [9]. 

Inflow region: m(£) <C m max , 7Ti'(£) > 0. In this region one can neglect the first 
two terms in (JJJ, as can be shown by a careful analysis of (|12|) . The remaining terms give 
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9U(m 2 )" — m'/m 2 . Integrating once, we find 9U(m 2 )' — —mT 1 + e" x by requiring the 
derivative m'(£) to vanish at the point with minimal slope e. Integrating again, we obtain an 
implicit equation for to(£), 

l8Ue [m 2 /2 + me + e 2 ln(m/e - 1)] = f . (10) 

Note that 1)10(1 is valid for to(£) > e + , to avoid the logarithmic singularity at to = e. In 
reality there is no singularity, because additional terms from (JJJ have to be included when 
TO — * £, which however are completely irrelevant in the remaining part of the inflow region. 

Scaling laws. - We are now prepared to investigate the scaling properties of large step 
bunches (M > 1). 

(A) The easiest is to find the size of the first terrace in the bunch l\, defined as in [11] by 
m{h = h ) = h /li, since this region is well described by l|10|) . Away from the singularity, for 
TO ^> e, i|10|) reduces to to(£) « ^/£j9Ue, which is of a similar form as the Pokrovsky-Talapov 
singularity found for static bunches [11]. This yields immediately 

h « (6Ue) 1/3 . (11) 

(B) To estimate the size of the minimal terrace in the bunch Z m ; n = l/m max , we multiply 
by m(m 2 y — 2m 2 m' — (2/3)(m 3 )' and integrate from £i to £ 2 : 



3U ,. 2 w\2 , fm 3 m 2 

— ((m 2 )') + b\ 

4 VV ' ' \ 3 2 



?2 

= -29. m 2 m' (£ + £ - h)d£+ - / - — '—d£. (12) 
J 3Jm 



6 

si €l 



Setting £i = 0, m(0) = e and £2 = £max with m(^ max ) = TO max , the left hand side gives 
ss fem^x/S for m max 3> 1. The second integral on the right hand side can be taken, noting 
that the function (m') 2 /m vanishes everywhere except in the narrow inflow region, where 

iffflf holds, yielding J^((m') 2 /m)d£ = [-to' - {iWm 2 )- 1 } ^ « (Sej/e 2 )^ 1 = 7 . To 
estimate the first integral, first set £2 = M in (|12[1 so that the left hand side vanishes due to 
the periodic boundary conditions. We obtain then Iq [0, M] = 2f2 J m 2 m' (£ + £0 — fe) d£ ~ 
Jo/3. Denoting by 7 = lim^/^oo Jo [0, £max] /In [0, M] the relative contribution to the integral 
from the segment [0, £ max ] for large M, we find from (fT2"|) dm^ ux 1=3 (1 — 7) Jq, or 



1/3 

1 ma X = ^min 



( 3QUe 2 b 
\ 1-7 



(13) 



Numerically we observe that the value of 7 indeed saturates to a fixed value for large M, and 
depends rather weakly on b and £/. Varying b and f around physically relevant choices of 
parameters, e.g., those in Fig^or in [11], we find 7 »s 0.7 ± 0.01. 

(C) Bunch width W. The definition of the bunch width depends on the convention used to 
assign steps to the bunch [11]. Here we define the bunch as the collection of terraces with sizes 
smaller than the mean terrace size I. We can obtain an estimate of W integrating the first term 
on the right hand side of (|12fl by parts: J Q M m 2 m! (£ + £0 — h) d£ = (1/3) J Q M m 3 (to — 1) d£, 
and arguing that the dominant contribution to the integral comes from the bunch interior. 
Then J^ 1 m 3 (to — 1) d£ ~ m 4 d£ — Qm^^W » M/ /(66), where Q < 1, and we have 
used f2 « 2>b/M. Substituting m max we get 
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For the parameters of Fig^ Q ~ 0.3237 for large Af , and it changes only slightly (±2%) 
under significant variations of 6 and U. 

(D) The minimal slope. As we have seen, many characteristics of the moving bunch are 
controlled by the single parameter e, which can be defined microscopically as the inverse size of 
the largest terrace in the outflow region between two consecutive bunches, ^ max = 1/e. In order 
to make a connection to earlier studies, we need the dependence of e on the bunch size M. 
From the asymptotic slope profile m w l/(36£) in the outflow region, one expects e w (36M) -1 
to leading order. Numerical studies suggest however strong finite size corrections even for large 
bunches, M < 200. The behaviour of e is rather well approximated by e _1 « 3baM — A(b, U) 
where a < 1 and typically A is not small, e.g., for b, U from FigQ] A w 13 in the range of 
bunch sizes 50 < M < 550. 

If we nevertheless use the asymptotic expression e w (36 Af ) _1 in ^ll [l -l)14 f we recover the 
scaling laws derived in [11] for static bunches, however with different numerical prefactors. 
This provides an a posteriori justification for the agreement between the predictions for static 
bunches and the numerical data in [11]. Noting that the dimcnsionlcss parameter S introduced 
in [11] is given by S = U/(bl) in the present units, we see that (fTTf) and (fT~5fl reduce to 
li ~ {2S/M) 1 / 3 and Z min w (13.3 x S/Af 2 ) 1 / 3 , which is to be compared with the expressions 
h = (4S/M) 1 ' 3 and l min = (16S/M 2 ) 1 / 3 for static bunches. From O we find that the 
ratio W/ (A/7 m ; n ) = [6Q(1 — 7)] -1 ~ 1.72, which is considerably larger than the corresponding 
number 1.29 in the static case. This reflects the fact that moving bunches are considerably 
broader than their static counterparts, because of the gradual increase of the terrace sizes in 
the outflow region, and explains the significant discrepancy between numerical and analytic 
estimates for W in [11]. Finally, we note that for a general step-step interaction exponent v 
in J3J) we arrive at generalized scaling laws l\ ~ , / min ~ £ 2 /i^+ 1 ) an d w ^ £ - 1 /(' y + 1 ) 

which are, apart from strong finite size corrections, consistent with the ones derived in [11]. 

Experimental considerations. - An important condition for the applicability of our con- 
tinuum theory is the smallness of the asymmetry parameter, 6 < 0.5. For step bunching 
induced by a conventional Ehrlich-Schwoebel (ES) effect during sublimation, 6 = (k + — 
/c_)/(fc+ + fc_), where k + (fe_) is the kinetic coefficient for attachment to a step from the 
lower (upper) terrace [11]; keeping 6 small then simply requires a weak ES effect. For current- 
induced step bunching in the attachment/detachment limited regime, one finds [14] 

_ kc cq a 2 FT e _ TFr e _ TgT e ml 

k^f ~ 2a 2 fc B T' ~ 2k B T 

where k = k + = denotes the attachment rate of adatoms to steps, a 2 the atomic area, 
c eq the equilibrium adatom concentration, F the electromigration force, r e the monolayer 
evaporation time, and g the step interaction strength. The quantity T = 2fcc cq a 4 is the 
mobility of an isolated step. The model Q of non-transparent steps is expected to apply in 
two of the four temperature regimes [4] in which step bunching is observed on Si(lll), around 
900° C and around 1250° C [16-18]. The material parameters given in [12,14] lead to the 
estimate 6 sa 14 in the low temperature regime and 6 ~ 0.3 in the high temperature regime. 
The latter is presumably an upper bound, since a recent estimate [20] of the kinetic length [7] 
Ik = D/k (where D is the adatom diffusion coefficient) at 1200° C indicates that the step 
mobility T increases less rapidly with increasing temperature than was assumed in [12]. The 
step interaction parameter U depends very sensitively on the mean miscut m = h Q /l. Taking 
the interaction strength to be g ~ 0.1 eV/A 2 at 1250° C, the parameters given in [12,14] yield 
U/l w 0.2 for I = 50 nm. 



(15) 
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Conclusions. - In this paper we have shown that the continuum equation JSJ faithfully 
represents the properties of moving step bunches for sufficiently small values of b, and that 
it can be used to extract accurate analytic expressions for the bunch shape and its various 
characteristic length scales. Our study reveals a central role of the point of minimal slope 
(or maximal terrace size) where the outflow region of one bunch joins the inflow region of 
the next. An important ingredient is a Fourier analysis of the discrete equations of step 
motion, which yields an inverse dependence of the bunch speed on bunch size under rather 
general conditions. Further work is needed to analytically derive the coefficient k in J5|, and 
to understand how to obtain the bunch speed directly from the continuum equation, possibly 
by exploiting a recently proposed analogy to front propagation problems [9] . In addition, it 
seems desirable to study the regime of large b, and to investigate the consequences of bunch 
motion for bunch coarsening beyond simple scaling arguments [9, 19]. 

The conditions assumed in this paper should be realizable in electromigration experiments 
on Si(lll) at temperatures around 1250° C. The gross features of the morphology, such as the 
maximal terrace size Z max and the bunch width W, could be followed in real time using REM, 
while for the more delicate measurements of Z m i n and h STM-studies would be preferable. 

* * * 
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